Tortuosity— porosity relation in the porous media flow 



00 

o 
o 

(N 
G 

(N 



>> 
i 

q 

O 

C/3 

Ph. 



> 

m 
rn 

o 

oo 
O 



Maciej Matyka, 1 Arzhang Khalili, 2,3 and Zbigniew Koza 1 

'institute of Theoretical Physics, University of Wroclaw, pi. M. Borna 9, 50-204 Wroclaw, Poland 
2 Max Planck Institute for Marine Microbiology, Celsiusstrasse 1, D-28359 Bremen, Germany 
Jacobs University Bremen, Campus Ring 1, D-28759 Bremen, Germany 
(Dated: February 2, 2008) 

We study numerically the tortuosity-porosity relation in a microscopic model of a porous medium 
arranged as a collectin of freely overlapping squares. It is demonstrated that the finite-size effects 
and the discretization errors, which were ignored in previous studies, may cause significant under- 
estimation of tortuosity. The simple tortuosity calculation method proposed here eliminates the 
need for using complicated, weighted averages. The numerical results presented here are in good 
agreement with an empirical relation between tortuosity (T) and porosity (0) given by T — I oc In <j>, 
that was found by others experimentally in granule packings and sediments. This relation can be 
also written as T — 1 oc RS/(j> with R and S denoting the hydraulic radius of granules and the 
specific surface area, respectively. 

PACS numbers: 47.56.+r,47.15.G-,91.60.Np 



I. INTRODUCTION 

In the low Reynolds number regime, flow through a 
porous matrix is governed by Darcy's law that links the 
fluid flux (discharge per unit area) q with the applied 
pressure gradient VP by the linear relation 



-VP, 



(1) 



where /i is the dynamic viscosity of the fluid and k is a 
proportionality constant known as permeability To 
a large extent, the proper description of the fluid flow 
through a porous medium depends on precise relations 
between the physical properties involved such as perme- 
ability and porosity (</>). In particular, much attention 
has been paid to deriving relations between k and tfi 0- 
In 1927 Kozeny developed a simple capillary model for a 
porous medium, and proposed the relation 



(2) 



where S is the specific surface area and Co is a dimension- 
less Kozeny constant that depends on the channel geom- 
etry 1]. Unfortunately. Kozeny's formula is not universal 
and does not hold for complicated porous geometries [l| . 
For example, it does not take into account pore connec- 
tivity and the fact that the specific surface area can be 
increased to an arbitrarily large value by removing some 
of the material to roughen the porous matrix surface in a 
fractal-like manner. On a purely physical ground namely, 
one would expect that removal of the material from a 
porous matrix would increase its permeability, whereas 
Kozeny's formula predicts just the opposite @]. 

One of the most widely accepted attempts to general- 
ize relation ([2]) was proposed by Carman [l], d, [f| , who 
noticed that the streamlines in a porous medium are far 
from being completely straight and parallel to each other. 



This effect can be described by a dimensionless parame- 
ter T called hydraulic tortuosity, 



T : 



(A) 
L ' 



(3) 



where (A) is the average length of the fluid paths and 
L is the geometrical length of the sample [l| . Using the 
tortuosity, Kozeny's relation ([3]) can be generalized to 



k = c 



T 2 S 2- 



(4) 



By fitting experimental data, Carman concluded that T 2 
is a constant factor (sa 5/2) over a wide range of porosi- 
ties. Later it was found that T 2 does vary with <f>, and 
can be as large as 50 for low porosity media [1, @] • 

Furthermore, it was realized that elongation of stream- 
lines not only affects the hydraulic discharge, but also me- 
diates other types of transport phenomena in the porous 
medium. This resulted in introducing several distinctive, 
experimentally measurable tortuosities obtained from a 
particular transport process, leading to diffusive [1, 
electrical 0, [13, [U and acoustic [ll| tortuosity defini- 
tions. There were also further theoretical attempts to 
define tortuosity [l], [f| [Hj]. However, all these tortuosi- 
ties, in general, differ from each other. Except for some 
very simple models d, H, [Hj], there is no clear consen- 
sus on the relation between these definitions. Among all 
these definitions, the one expressed in Eq. ([3]) is not only 
the simplest, but also widely adopted in theoretical stud- 
ies, for it ties tortuosity with the underlying geometry 
and topology of the porous medium. 

It has been known since long that flow through a 
porous medium depends on many factors such as poros- 
ity, tortuosity, granule shape and size distribution, satu- 
ration, Reynolds number, etc. For proper understanding 
of transport phenomena in porous media, however, it is 
essential to depart from simple systems with a limited 
number of well-defined control parameters. Therefore, 



in this paper we investigate the hydraulic tortuosity (as 
defined in Eq. ((SJ) in a creeping flow through a porous 
region constructed by a two-dimensional lattice system 
with a uniform, randomly distributed and freely overlap- 
pin g so lid squares. This model, first used by Koponen et 
al. [1 21 ] . is simple enough to allow a numerical solution 
with the advantage of having porosity (0) as the only 
control parameter. It is particularly suitable for study- 
ing tortuosity-porosity relation, as at high porosities the 
streamlines are almost straight, whereas for low porosi- 
ties they become wiggly. 

Using the Lattice Gas Automata (LGA) method, Ko- 
ponen et al. (l2| solved the flow equations for a porosity 
range of <p £ [0-5, 1], and concluded that 



T = p(l 



1, 



(5) 



where p is a fitting parameter. However, later they found 
the relation not being consistent with the results obtained 
for the porosity range d> & [0.4,0.5], and suggested [HI 
to replace relation §5§ by 



T = p 



(1 



(6) 



in which <p c w 0.33 is the percolation threshold while p 
and m are some empirical parameters. Still, as an ad hoc 
formula with two fitting parameters, this relation can not 
be considered a universal law. In addition, the data used 
to derive ([6]) suffer from systematic errors, as neither the 
impact of a finite system size nor the effect of the space 
discretization were taken into account. 

The aim of this paper is to cary out a numerical simu- 
lation for analysing the tortuosity — porosity relation in a 
system of freely overlapping rectangles with a high accu- 
racy. In addition, we provide a simplified algorithm for 
T calculation without the need for implementing compli- 
cated, weighted averages of streamline lengths. 

The structure of the paper is as follows. Section |TT] 
specifies the model and the numerical techniques used. 
Special attention is paid to the description of the non- 
trivial numerical technique for the tortuosity. Main re- 
sults, including a detailed numerical error and finite-size 
analysis are provided in Sec. IIHI Finally, the results are 
discussed in Sec. IIVI 




FIG. 1: An example of two 800 x 800 porous matrices con- 
structed by randomly placed and freely overlapping rectangles 
of size 10 x 10 for two different porosities (/>■ 



aligned with the y-axis of an x-y Cartesian coordinate 
system to model the gravity. Following previous works 
[2 EH i periodic boundary conditions have been imposed 
in both directions. 

Two examples of such porous systems are depicted in 
Fig. [TJ The dark areas represent fixed solid obstacles, 
while the white part is occupied by the fluid. 



B. Numerical techniques 

Numerical solution of the model defined above con- 
sists of five main steps: (i) generation of a porous matrix 
of a known porosity; (ii) solving the flow equations in 
the low Reynolds number regime; (iii) finding the flow 
streamlines; (iv) determining the tortuosity of the flow; 
(v) error analysis. 



1. Construction of the porous matrix 

A porous matrix of a given porosity <p can be generated 
using the method of (l2l . [T3 |. Starting from an empty 
system, solid squares are added at random positions until 
the desired porosity has been reached. The porosity is 
calculated as the fraction of unoccupied lattice nodes. 



II. MODEL 

A. General description 

The system of interest consists of a square lattice 
(L x L) in which a number of solid squares (a x a lat- 
tice units) have been placed at random locations to form 
a porous matrix (1 < a <C L). The squares are fixed 
in space but free to overlap. The only restriction is that 
their sides must coincide with the underlying lattice. The 
remaining, void space is filled with a liquid. The con- 
stant, external force imposed on the porous medium is 



2. Lattice Boltzmann method for solving flow equations 

To solve the flow equations, we applied the Lattice 
Boltzmann Model (LBM) (l5| with a single relaxation 
time collision operator [161 ]. This method proved useful 
in microscopic model simulations of flow thro ugh p orous 
media for various conditions and flow regimes 2, 17]. It is 
a numerical technique that rests on the Boltzmann trans- 
port equation discretized both in time and space, and is 
expressed in terms of the velocity distribution functions 
rii in the form 

n t + 1 {r + c, i ) = n\{r) + A\{r)+F l , (7) 
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FIG. 2: A schematic view of the system before (a) and after 
(b) grid refinement. In this example the system size is 4 x 6, 
the porous matrix is made of three quads of size 2x2 (A, B 
and D), and the refinement level fc ro f = 2. 



where 



0.. 



.8 identifies lattice vectors c^, t is an 



integer time step, r denotes a lattice node, A|(r) is the 
collision operator at r, and Fj represents the i-th compo- 
nent of the external force. We used the time unit equal 
to the relaxation time, which yields the kinematic vis- 
cosity v = 1/6 [15]. This, in turn, simplifies the form 



of the collision operator A*(r) 



1 (r) — n\(r) with 



n° q (r) being the equilibrium value of rii at r. The ex- 
ternal force was taken into account using a method of 
Ref . [HI : half of the momentum was transferred directly 
into the equilibrium distribution function during the col- 
lision step, whereas the other half was included into the 
transport equation. Because we were interested in the 
solution of a slow, laminar flow, we utilized the equilib- 
rium distribution function n- q linearized in the velocity 
as 



n° q = Wip [1 + 3(it • Ci)] , 



(8) 



in which u is the macroscopic velocity vector and Wi are 
some weighting coefficients that depend on the lattice 
structure and dimension [lij [2(| . 

One problem with the LBM method is that it is inca- 
pable of resolving the macroscopic Navier-Stokes equa- 
tions for channels narrower than about 4 lattice units 
(l~5| . This limitation becomes particularly important at 
low porosities, for which the number of very narrow pas- 
sages increase enormously. To bypass this problem, a 
standard numerical mesh refinement procedure was used. 
Starting from the original lattice taken to generate the 
porous matrix, each of its L 2 elementary quads were sub- 
divided into k rc { x fc rc f smaller quads with fc rc f = 1,2,... 
being the refinement level. The resulting k TC fL x fc re f£ 
computational grid of vectors r in Eq. ([7]) will then be 
formed from the centers of the small quads. With this 
choice, the identification of the interface between the 
porous matrix and the free space is facilitated. 

A schematic picture of the refinement procedure is 
shown in Fig. [5] Note that the refinement effectively in- 
creases the number of the lattice nodes between any two 
points by the factor fc re f, and that the smallest channel 
width is fc re f + 1 lattice units. 

After initialization, the LBM computational loop of 
advection and collision continued for 5 x 10 3 time steps. 
By using the mid-grid bounce-back rule applied to the 




FIG. 3: Velocity magnitudes squared (u 2 = u 2 , + v 2 ,) calcu- 
lated on a 600 x 600 numerical grids, which corresponds to a 
200 x 200 lattice (refinement level fc ro f = 3). The square block 
sizes were 10 x 10 (i.e., 30 x 30 after refinement) and the poros- 
ity was cj> — 0.65. Periodic boundary conditions were assumed 
in both directions. The grey squares in the figure represent 
the solid part of the medium, whereas the black region is the 
porespace open to fluid flow. 



no-slip boundaries, second order accurate solutions, both 
in space and time, could be achieved (l5| . An example 
of the velocity field calculated with this method for a 
low-porosity medium is shown in Fig. [3l 



3. Flow streamlines 

After obtaining the velocity field u at each grid point, 
streamlines could be obtained by solving the equation of 
motion for the trajectories r(t) of massless particles [2~lj |. 



dr 



u(r). 



(9) 



The u(r) values of points lying between two grid nodes 
were obtained using bilinear interpolation. Due to com- 
plex boundary conditions and extreme velocity differ- 
ences on the grid, the 4-th order Runge-Kutta algorithm 
with adaptive time stepping was used [22l ]. 



4- Tortousity 

The tortuosity is defined by Eq. as the ratio of the 
average length of all particle path lines passing through a 
given cross-section during a unit time period to the width 



of the sample [![ leading to 

T= 1 I A u y( x )H x )dx ^ 
L I A u v( x ) dx 

in which A is an arbitrary cross-section of the system 
parallel to the x axis, X(x) is the length of the streamline 
cutting A at x, and u y {x) is the y component of the trial 
particle velocity at x. 

The integrals in (TIT))) have been obtained in the liter- 
ature either by the Monte Carlo integration [12, EH H3| 
or by direct quadratures @. In the former method, the 
lengths of the streamlines passing through randomly cho- 
sen points within the pore volume are averaged using 
proper weights. In the latter method T is approximated 
by the relation 

T w i T,jM x j)Hxj)^x 3 

~ L V, U y {x j )Ax j 

where Axj = Xj+\ — Xj are discretization intervals of 
A. In principle, both approaches should yield the same 
results, but both can be easily misused. For exam- 
ple, some researchers used the Monte Carlo integration 
with streamlines passing through points chosen randomly 
from a uniform distribution over the whole pore space 
some others calculated streamlines cutting all 
lattice nodes [HI, whereas others recorded all stream- 
lines crossing every lattice node along a chosen inlet plane 
0, [l3j]. However, such 'uniform' approaches are not co- 
herent with the reality of low porosity systems, in which 
transport is mostly carried out only through few 'con- 
ducting' channels (cf. Fig. [3]). Consequently, the sums 
in (flTj) contain most probably many terms of practically 
negligible magnitudes. To avoid this problem, we used 
Eq. (|11[) with a constant-flux constraint between two 
neighboring streamlines, 

u y (xj)Axj = const. (12) 

With this choice, Eq. (jlip immediately simplifies to 

1 1 N 

T «z^E A ^-)' ( 13 ) 

where N is the number of the streamlines generated. 
Note that all terms in this sum are of the same order 
of magnitude. 

Thus, to calculate T, a horizontal cross-section A is 
chosen. Next, the coordinates of the initial points x j are 
determined using (|12[) , and the corresponding streamlines 
are found by solving ([9]) in both directions until the so- 
lutions hit the system edges (y — and y = L). Finally 
their lengths are plugged into (p~5]) . 

It should be noted that not all streamlines passing 
through Xj cut both horizontal edges y — and y = L 
(see Fig. This may happen if a streamline passes 
through a region with extremely low fluid velocity. There 
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FIG. 4: Streamlines generated with the constant-flux con- 
straint (|12|) for the same system as in Fig. [3] (N = 59). The 
horizontal line represents the cross-section y = L/2 on which 
the initial points Xj were chosen. Only those streamlines are 
shown for which \(xj) could be determined numerically. The 
arrows point at regions discussed in the text. 



are two types of such regions: dead-end pores and vol- 
umes where the fluid stream is split by an obstacle (or 
merges behind it). An example of Xj located in a dead- 
end pore is shown by arrow (a) in Fig. |4j Arrow (b) 
points at Xj corresponding to an incomplete streamline 
that flows from a region (c) where the stream merges 
after being split by an obstacle. 

To bypass this problem, in calculation of the sum (fT3"|) 
only those streamlines were taken into account, whose 
lengths Xj can be determined. The error induced by this 
procedure is discussed in Sec. IIIII 



5. Error analysis 

Tortuosity values (T) calculated directly from (fTB")) con- 
tain errors arising from different sources. While statisti- 
cal errors result from randomness in the porous matrix, 
discretization errors appear when approximating the in- 
tegrals in (flT)]) by finite sums, and when solving flow equa- 
tions by discrete lattices. Finite-size errors could emerge 
also as a consequence of approximating a macroscopic 
system with a microscopic model. Details of the error 
analysis are addressed in the next section. 
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FIG. 5: The basic quantities that affect tortuosity calculation 
by Eq. (|10p : (a) velocity component u y (x); (b) local tortuos- 
ity t(x) — \(x)/L; (c) the product u y (x)r(x); (d) the ratio 
of the minimum to maximum speeds along streamlines. All 
quantities were determined for the system shown in Fig. [4] 
with the cross-section y — L/2 and the uniform discretization 
Axj = L/N (N = 1200). The tortuosity for this system is 
T ~ 1.45. 



III. RESULTS 

To begin te discussion, the structure of the integrands 
in (flO]) is examined. Figure [5] shows u y (x), the local tor- 
tuosity t(x) = \(x)/L, their product u v (x)t(x), and the 
ratio of the minimum to maximum trial particle speeds 
for the same system as in Figs. [3] and HI The data for 
this figure were obtained from the streamlines originat- 
ing at the cross-section y — L/2. As expected, the ve- 
locity profile is continuous and piecewise-differentiable, 
and partially resembling that of a Poiseuille (parabolic) 
flow. The negative value of u y near x = 160 indicates 
that some streamlines are cutting the initial cross-section 
many times. This effect must be taken into account to 
avoid multiple counts of the same streamline in Eq. (|13p . 

In contrast to u y , the local tortuosity r is a discontin- 
uous function of x (Fig. [SJa). Each jump in the t(x) plot 
corresponds to the stream flow splitting into (or merg- 
ing from) two parts upon meeting an obstacle. Thus, 
for a finite-size system, t(x) is a piece-wise continuous 
function with a certain number of discontinuities, equal 
to the number of ,, islands" existing in the porous do- 
main. Consequently, the product u v (x)t(x) is also dis- 
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FIG. 6: Tortuosity T as a function of the system size L for 
a = 10 and <j> — 0.7 and 0.5, averaged over M = 30 samples. 
The solid line is the best fit calculated with Eq. The 
error bars represent the standard error of the mean. 

continuous (Fig. [5t). Moreover, the problem of finding 
the coordinates of discontinuity points is numerically ill- 
conditioned. These two factors greatly complicate the de- 
termination of the enumerator in Eq. (|10p. and introduce 



an additional source of errors in (1131). For i, near a dis 



continuity point, even small numerical errors may result 
in a significant jump in X(xj). Two countermeasures were 
taken to reduce the impact of this phenomenon, which is 
closely related to the problem of "missing streamlines" 
discussed in Sec. Ill Bl First, a check is made to find 
out how T calculated from (flU)) varies with N. Here, an 
optimal value of N w L was found. Second, the tortuos- 
ity was always calculated as an average over 8 different 
cross-sections. This approach not only reduced the error 
resulting from approximating (JTUJ) by (flU)) . but also gave 
some estimation on its magnitude. The errors were found 
to be maximum for low porosities, but even for <p = 0.45, 
the relative error was less than 0.5%. 

The large number of discontinuities in t[x) implies 
that the fluid velocity along a typical streamline may 
vary by many orders of magnitude. This is shown in 
Fig. [5]i, in which the ratio of the minimum to maximum 
fluid speed (w m i n /u max ) along the streamline cutting the 
cross-section y = L/2 at x is plotted. In this particular 
case Umin/nmax < 0.26 and drops to at all positions 
where t(x) became discontinuous. A comparison of pan- 
els (a) and (d) in Figure [5] shows that streamlines passing 
through a region with relatively high fluid velocity will 
likely hit regions where the fluid is almost stagnant. For 
this reason it is essential that Eq. be solved with a nu- 
merical method that uses variable step lengths and local 
error control. 

Next, finite-size effects were analysed. Figure [S] shows 
the dependency of the tortuosity T on the system size 
L for two porosities <f) — Q.7 and <fi — 0.5 averaged over 
M = 30 samples. The lines represent fits to 



T(L) = Too — feexp(-cL), 



(14) 
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FIG. 7: The dependency of T on the system size L for (f> = 
0.45 and three refinement levels fc rc f (symbols). The lines 
represent best fits to the function given by Eq. (I14|l . 



FIG. 8: The dependency of the tortuosity T on the porosity 
4>. Our data obtained with the LBM method and Eq. (|13[) 
(circles); relation (|SJ derived for the same system by Koppo- 
nen et al. (dashed line); the best fit to Eq. ()15b|l (solid line). 



where T^, b and c are free parameters. These, as well 
as many other fits, not shown here, proved to be suit- 
able enough to estimate the tortuosity of an infinite sys- 
tem Tqq. This procedure also enabled us to estimate a 
characteristic system length L* above which T does not 
change significantly with L. It turned out that L* « 200 
for (j) > 0.55 and L* w 300 for 0.45 < (j) < 0.55. In all 
cases analysed, T was found to be an increasing function 
of L. Thus, it becomes clear that ignoring finite-size ef- 
fects and using L < L* would lead to an underestimation 
of T. 

Following this, the sensitivity of tortuosity measure- 
ments to the numerical mesh refinement was examined. 
Figure [7] depicts the values of T(L) for three refinement 
levels k rc f = 1,3,4. We found that significantly de- 
pends on fc re f only for fc re f < 3. This threshold value is in 
accord with the criteria, mentioned above, for the LBM 
method to reconstruct the Navier-Stokes equations [l5| . 
Note that ignoring numerical mesh refinement would lead 
to a significant underestimation of T. This shows that 
narrow passages are a relevant factor that affects trans- 
port properties of a porous medium. Numerical mesh 
refinement is thus particularly important at low porosi- 
ties, where narrow passages are common. 

After finding the minimal requirements on the mesh re- 
finement level fc rc f and the system size L* , the tortuosity- 
porosity relation could be determined. For a given 
4> = 0.45, 0.5, . . . , 0.95, a system size of L = L* and a 
refinement level of k rc { — 3 were chosen. At each <f>, T 
was calculated for M porous matrices, with M ranging 
from 30 (for = 0.95) to 140 (for = 0.45), and the 
results are shown as circle symbols in Fig. [51 For com- 
parison, we also plotted the best-fit curves obtained for 
exactly the same system by Koponen et al. [ijj], see Eq. 
([6]). Obviously, the results of Koponen et al. lie signif- 
icantly below those obtained by us. This is due to the 
fact that in the work of Koponen et al., a rather small 
system (L < L*) is considered without numerical mesh 



refinements. Hence, the discrepancy can be explained 
as a consequence of finite-size effects and discretization 
error analysis, which are missing in their study. 

We fited our data to four tortuosity-porosity relations 
proposed by other researchers: 



T(4>) 
T(<f>) 

T(<P) 



1 — p In (f>, 
l+p(l-0), 



(15a) 
(15b) 
(15c) 
(15d) 



where p is a parameter. The first of them was proposed 
for the electric tortuosity by Archie (1942) [24[. The sec- 
ond equation (with p = 1/2) was found in a theoretical 
study on diffusivity of a model porous system composed 
of freely overlapping spheres by Weissberg (1966) [25| . 
The same relation (with p w 0.86 and p s=s 1.66) was 
also reported in measurements of the hydraulic tortuos- 
ity for fixed beds of parallelepipedal particles with differ- 
ent thickness-to-side ratios by Comiti and Renaud (1989) 
(26| and in recent measurements of electrical tortuosity in 
fixed beds and suspensions of glass spheres by Barrande 
et al (2007) [27[. Equation (|15c[) is an empirical relation 
found for sandy (p — 2) or clay-silt (p = 3) sediments 
by Iversen and J0rgensen (1993) [28] . Finally, equation 
(I15d[) . with p = 32/97T w 1.1, was recently obtained in 
a model of the diffusive tortuosity in marine muds by 
Boudreau and Meysman (2006) . 

Even though we treated p in all these formulas as an 
adjustable parameter, only Eq. (|15b[) gave a satisfactory 
fit for p = 0.80 ± 0.01 (the reduced chi-square statistics 
« 0.8), a value comparable with those of Comiti and 
Renaud HI (p ss 0.86 and p w 1.66). This fit is plotted 
in Fig. [5] as a solid line. 

The fact that our system obeys Eq. (|15b[) has a rather 
interesting and unexpected consequence. As shown pre- 
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viously [TH, the specific surface area S in the model of 
freely overlapping squares of side a satisfies the relation 

S=~^hx4>, (16) 

with R — a/2 denoting hydraulic radius of obstacles. 
With this, equation (|15bj) simplifies to 

T-lociZ?. (17) 

<P 

IV. DISCUSSION AND CONCLUSIONS 

As shown by the present study, obtaining hydraulic 
tortuosity from numerical simulations contains many hid- 
den problems, which may lead to incorrect conclusions. 
When a fluid stream hits an obstacle, it splits and then 
merges, causing a discontinuity in streamlines. The 
bounding streamline of each obstacle separates the two 
splitting (or merging) streams. The location of such 
streamlines is a priori not known, and the problem of 
finding the streamlines within those regions is numeri- 
cally ill-conditioned. If the system is sufficiently large, it 
is inevitable that majority of streamlines pass through 
such "ill-conditioned" regions. Moreover, the velocity 
magnitude along the streamlines can vary by many orders 
of magnitude. 

These problems resemble those encountered in an- 
other type of computer simulations, molecular dynamics, 
where discontinuities (as well as large variations in veloc- 
ity) are caused by collisions. After several such events, 
the computer-generated particle trajectories have liter- 
ally nothing to do with the exact solutions. Still, molec- 
ular dynamics is one of the most successful methods of 
computer physics. The reason is that physical quanti- 
ties never depend on exact trajectories of the individual 
particles — it is sufficient to ensure that the solution keeps 
on a constant energy surface. This analogy makes us be- 
lieve that despite all difficulties, hydraulic tortuosity is 
a well-defined quantity that can be reliably obtained by 
numerical methods. All simulations we performed with 
different cross-sections, different numbers of streamlines, 
different choices of streamline starting points, different 
numerical ordinary differential equation solvers, resulted 
in almost the same numerical values. This implies that — 
just as in molecular dynamics — small local errors, which 
are unavoidable in computer simulations, are of marginal 
importance. This corresponds to the fact that in real 



fluid flow, trajectories of individual molecules are not lim- 
ited to single "theoretical" streamlines, but are affected 
by diffusion at low velocities or turbulence at higher ve- 
locities. 

It is also interesting to mention that the problem of 
finding T is easiest at high porosities, where nearly each 
obstacle constitutes a separate island. Although the 
number of discontinuities is very large, in general they 
tend to average out. At low porosities, however, se- 
vere problems may arise as discontinuities are much fewer 
in number (which means no "averaging out") but much 
larger in magnitude (which results from increased island 
sizes). With this study we have demonstrated how sen- 
sitive tortuosity computations are to finite-size effects, 
discretization errors and large variation of fluid speed 
along streamlines. The system size must be large enough 
to ensure development of chaotic "splitting and merg- 
ing" flows, that are characteristic of real granular sys- 
tems. Discretization errors creep into the system from 
several places, most notably in narrow channels, and can 
be avoided by numerical mesh refinement. Our results 
concerning large fluid velocity variations along stream- 
lines are a clear indication for revising those tortuosity 
definitions which assume a constant fluid velocity along 
a streamline @, [l3| . They also show that numerical de- 
termination of streamlines requires using advanced nu- 
merical integrators with adaptive step lengths and local 
error control. 

When streamlines are generated using the constant- 
flux constraint, the tortuosity can be calculated simply 
as an average over the streamline lengths. This method 
reduces the computation errors and does away the need 
for determining the streamlines in dead-end pores. 

The numerical data presented in this study were 
found to be in good agreement with those generated by 
Eq. (|15b[) . obtained from previous experimental studies. 
However, this relation cannot be applied close to the per- 
colation threshold, where tortuosity diverges. Note also 
thatEq. (fi~5bf was found to describe both hydraulic and 
diffusive tortuosities, i.e. quantities that are certainly cor- 
related, but in a way that has not been well established 
yet. It is not clear whether our findings reflect equiva- 
lence of these quantities, or whether they can be consid- 
ered as a coincidence. 

We also found that in the model of freely overlapping 
squares, a very simple relation (jTTJ) holds between tortu- 
osity, porosity and the specific surface area. This equa- 
tion is closely related to Eq. (|15b[) , and it expected to be 
valid for all systems, for which Eq. (|15b[) holds. 
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